arXiv:1508.05409vl [astro-ph.IM] 21 Aug 2015 


A New Diagnostic of Magnetic Field Strengths in 
Radiatively-Cooled Shocks 

P. Hartigan A. Wright 1 

ABSTRACT 

We show that it is possible to measure Alfvenic Mach numbers, defined as 
the shock velocity in the flow divided by the Alfven velocity, for low-velocity 
(V shock 100 kms -1 ) radiative shocks. The method combines observations of 
bright forbidden lines with a measure of the size of the cooling zone, the latter 
typically obtained from spatial separation between the Balmer emission lines 
and the forbidden lines. Because magnetic fields become compressed as gas in 
the postshock region cools, even relatively weak preshock magnetic fields can 
be detected with this method. We derive analytical formulae that explain how 
the spatial separations relate to emission-line ratios, and compute a large grid 
of radiatively-cooled shock models to develop diagnostic diagrams that can be 
used to derive Alfvenic Mach numbers in flows. Applying the method to existing 
data for a bright knot in the HH 111 jet, we obtain a relatively low Alfvenic 
Mach number of ~ 2, indicative of a magnetized jet that has super-magnetosonic 
velocity perturbations within it. 


Introduction 


Supersonic flows are common throughout the interstellar medium and govern the dy¬ 
namics of a variety of objects such as stellar jets (Frank et al. 2014) and supernovae remnants 


e.g. Reynolds 2008), and also play an important role in affecting the observed spectra in 


protoplanetary nebulae ( |Riera et al. 2003) and novae shells ( Williams||2013 ). In cases where 
shocked gas cools by radiating emission lines, line ratios provide measurements of the densi¬ 


ties, temperatures, ionization fractions, and abundances throughout the flow (e.g. Lazendic 


et al. 2006). Emission line widths are a powerful way to quantify the amount of turbulence 


and nonthermal motions present in these objects, and high-resolution images that detect 
proper motion have begun to open up the time-domain, making it possible to observe shear, 
turbulent mixing, fluid-dynamical instabilities and cooling in real time (Hartigan et al. 2011). 
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This wealth of information not only defines the physical characteristics at each location 
in flows and helps to clarify the nature of the driving sources, but can also identify the most 
important atomic or plasma processes that operate in supersonic flows. For example, the 
presence of both a narrow and a broad velocity component in the Hct emission at the blast 
wave of many supernovae remnants arises from the high cross section of charge exchange 
between postshock ions and neutrals (Chevalier et al. |l9~80). The resulting Ha emission-line 
shapes constrain the velocity and orientation of the shock front, and the degree to which ion 
and electron temperatures are equilibrated behind the shock (Ghavamian et al. 2001). The 


relatively large line width associated with the narrow component points to the possibility 
of a cosmic-ray precursor (Raymond et al. 2011). In the case of stellar jets, H II regions 
visible ahead of the strongest shocks define the regime where radiative precursors begin to 
affect shock dynamics, and the presence of transient bright spots where shock waves intersect 
leads to the physics of Mach stem formation and evolution (Hartigan et al.||2011 Frank et al. 


2014). 


Although magnetic fields can play an important role in shocked flows on both large and 
small scales, especially for systems like stellar jets where the driving source is likely to be 
an MHD disk wind of some sort, it is difficult to infer field strengths in stellar jets directly 
from optical emission line spectra. Optical lines show no discernible Zeeman splitting, and 
emission line ratios in shocked gas are degenerate with respect to field strength in the sense 
that spectra from low-velocity shocks without a preshock magnetic field closely resemble 


those from a higher-density, higher-shock velocity object with a field (Table 7 of Hartigan 


et al. 1994). If a bow shock is well-resolved spatially and emits in [O III] near its apex, 
one can use the extent of the [O III] emission and the shape of the bow to infer the shock 
velocity, and the Ha flux to determine the preshock density. With this information in hand, 
the observed electron densities from an emission line ratio like [S II] A6716/A6731 give the 
compression, and hence the strength of the preshock and the postshock field that is oriented 


in the plane of the shock (i.e., perpendicular to the normal of the shock front; Morse et al. 


1992, 1993) 


Most jets do not emit detectable synchrotron emission, but nonthermal radio continua 
have been measured in a few bright knots, which are typically unresolved spatially. Fields 
have also been inferred from synchrotron continuum flux levels and from polarization mea- 


]Ray et al. 

1997 

Chrysostomou et al. 2007; 

Carrasco-Gonzalez et al. 


2010). However, converting these measurements to an estimate of the field strength requires 


knowledge of the electron density in the radio-emitting region immediately behind the shock 
where the density is changing rapidly, and also involves assumptions about the filling factor 
of the gas, the degree of ionization, and the power spectrum of the electrons. 
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Magnetic field strengths inferred from the resolved-bow-shock method and from con¬ 


tinuum polarization both indicate surprisingly weak fields. Morse et ah (1992, 1993) found 
magnetic fields of ~ 30/rG in the preshock gas ahead of bright bow shocks in HH 34 and 
HH 111. The observed compression ratio between the preshock gas and the radiating gas is 
~ 40 (a factor of ~ 4 at the shock and another ~ 10 in the cooling zone), so that the field 
strength is ~ lmG where forbidden lines emit at temperatures of ~ 7000K. The electron 
densities and ionization fractions determined at that location in the postshock gas imply a 
total density of ~ 3 x 10 4 cm -3 , resulting in an Alfven speed of only ~ 10 kms" 1 in the 
dense cooling zone. Hence, while the magnetic signal speed is comparable to the thermal 
velocity within shocked knots in the flow, it is negligible compared with the bulk flow ve¬ 
locities of ~ 300 kms -1 . Field strengths inferred from radio continuum observations of the 


unusually strong shocks in HH 81 are also only 100/iG (Chrysostomou et al. 2007). These 


results are in agreement with proper motions of jets like HH 111 (oriented in the plane of 
the sky), where internal velocity variations on the order of 40 kms -1 produce shocks and 
not magnetic waves, implying that the magnetic signal speeds in the jet are at most a few 
tens of kms -1 (Hartigan et al. 2001). 


Two factors help to account for the anomalously low fields observed in stellar jets. 
First, a pulsed flow with supersonic velocity perturbations will collect magnetic held into 
knots, while rarefactions in the how stretch the held out. This effect lowers the magnetic 
signal speed ahead of jet knots and allows smaller velocity perturbations to form shocks that 


otherwise would not do so (Hartigan et al. 2007). Another way to reduce the held in the jet 


is to remove it by reconnection. Some evidence exists for reheating at tens of AU from the 
sources both from high-spatial resolution observations of optical emission lines made with 


HST (Hartigan & Morse 2007), and X-ray observations of quasi-stationary knots that are 


offset spatially from their sources (Schneider et al. 2011). In laboratory settings, magnetic 


towers become kink-unstable close to their driving sources ( 

Moser & Bellan 

2012 

), and the 

same thing happens in some numerical simulations of disk winds ( 

Staff et al. 

2015 

). Kink 


instabilities could result in reconnection as the held geometry becomes increasingly complex. 

With the importance of measuring held strengths in MHD hows self-evident, it is worth¬ 
while to revisit the subject to see whether it may be possible to improve upon optical 
emission-line analysis of jets. Optical emission lines make the best connection to the source 
of any MHD wind, because they arise within the material driven directly from the source. 
In contrast, molecular emission lines in these outhows are typically slower and wider, and 
often come from a sheath of material entrained by the jet. It has been known for almost two 
decades that Ha and [S II] emit from different locations in jet shocks, with Ha coming from 
collisional excitation immediately behind the shock and [S II] following in a cooling zone 
(Heathcote et al. 1996). Surprisingly, this observational constraint has not been used to help 
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infer field strengths, despite the fact that even a weak preshock magnetic held affects this 
spatial offset in a significant way (e.g. compare models B100 and E100 in Hartigan et ah 


1987). 


In this work we show that the additional constraint of a spatially-resolved cooling-zone 
breaks the degeneracy of the line ratio information and makes it possible to estimate Alfvenic 
Mach numbers of radiative shock waves stellar jets directly, using only quantities that are 


readily observable. We discuss the physics of magnetized cooling zones in Sec. |2.1[ and present 
a large grid of shock models and their resulting diagnostics in Sec. PI Although the ideal 
observations have yet to be made, we collect the best available data to obtain preliminary 
estimates of the Alfvenic Mach number for a shock in the HH 111 jet in Sec. [3j and conclude 
with a short Summary section. 


2. Measuring Alfvenic Mach Numbers Using Emission Line Ratios and 

Cooling Distances 

2.1. Physics of Magnetized Cooling Zones 

As gas enters a shock front, the pressure, temperature, and density immediately behind 
the shock are determined by the Rankine-Hugoniot equations, derived by balancing the mass, 


momentum, and energy fluxes on either side of the shock (e.g. Draine 1993). Material that 


has passed through the shock is heated as the bulk kinetic energy in the flow transforms 
into thermal energy. The hot gas then cools by emitting line- and continuum-radiation, and 
within this cooling zone the temperature decreases and the density increases in response to 
the energy lost by radiation. The energy equation in the cooling zone of a partially-ionized 
steady-state shock is easily described by a formalism that relates the enthalpy per nucleus in 
the gas to the radiative energy loss rate (Raymond 1979"). Fig. [l] depicts an example of the 
spatial distribution of emission-lines in the cooling zone of a typical radiative shock front. 

Several types of magnetized shocks are possible depending on the orientation of the 
held relative to the shock front (e.g. Chapter 7 of Gurnett & Bhattacharjee 2005). In 


general, one must solve a cubic equation known as the shock adiabatic to obtain the jump 
conditions for density and pressure. If the held is oriented perpendicular to the shock normal 
(perpendicular to the velocity vector and parallel to the plane of the shock), there is a single 
real solution. In this case, an MHD shock forms as long as the incident velocity exceeds that 
of the fast magnetosonic speed V/ ms = [V^ + C|] 1/,a , where C 5 is the sound speed in the 
preshock gas and Y a is the Alfven speed in the preshock gas. If the direction of the held 
is aligned parallel to the normal to the plane of the shock (i.e., along the velocity vector) 
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two solutions are possible, one where the jump conditions are the same as those of the B=0 
hydrodynamical case, and another so-called switch-on shock where surface currents impart 
a tangential component to the field in the postshock gas. 

For the general case where the field is at an angle to the surface, one can transform 
into the deHoffmann-Teller frame where v || B, and the equations have slow, intermediate, 
and fast solutions that correspond as the compression approaches unity to, respectively, the 
slow MHD wave speed, the Alfven speed, and the fast MHD speed. Once the Alfvenic Mach 
number of the flow > 2, only the fast solution remains. Throughout this paper we have 
used the fast MHD solution, and when we refer to a magnetic field strength we mean the 
component perpendicular to the flow that lies in the plane of the shock. This component 
should be the dominant one for a typical jet geometry because in an MHD disk wind at 
locations far from the source, fields in jets should be mostly toroidal and therefore lie in 
the plane of the shock for a velocity-variable flow unless turbulence or kink instabilities 
randomize the direction. 

The addition of magnetic fields alters the structure of the cooling zone in important 
ways. Because even small ionization fractions (> 1CT 4 ) suffice to tie the magnetic field to 
the gas, the magnetic field strength B is proportional to the gas density n throughout the 
entire cooling zone. The total pressure in the cooling zone is roughly constant, so n ~ T~ 4 
as long as thermal motions dominate the pressure. However, the magnetic pressure B 2 /87r 
~ T~ 2 , so as the temperature in the cooling zone drops, the ratio of the magnetic pressure 
to the thermal pressure rises, and typically exceeds unity in the forbidden-line zone (see 
Figs. [2j [ 3 ]). Thus, even if the preshock magnetic pressure is orders of magnitude below the 
ram pressure at the shock front, by the time the flow radiates forbidden lines, the field 
usually dominates the pressure. A strong field in the cooling zone affects the observed line 
radiation mainly by (i) lowering the density, and (ii) by expanding the size of the cooling 
zone. 


From an observational standpoint, one can easily measure the electron density in the 
cooling zone from, for example, the red [S II] ratio A6716/A6731, but this measurement alone 
does not constrain the field strength because the density in the cooling zone also depends 
linearly upon the preshock density. The shock velocity also matters because higher shock 
velocities produce hotter postshock temperatures, which result in a larger temperature drop 
and subsequent larger density increase once the gas reaches the forbidden line zone. In 
addition, the preshock ionization fraction of the gas also has a small effect on cooling zone 
densities because more neutral preshock conditions result in lower postshock temperatures 
as some of the ram energy goes into ionizing the postshock gas (Cox & Raymond 1985). In 
practice, many emission lines occur in the cooling zone, and combining all of them narrows 
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the possible parameters considerably. However, as noted in the Introduction, one cannot 
infer magnetic field strengths from emission line ratios alone because adding a preshock field 
affects emission line ratios in a similar manner to lowering the preshock density and shock 
velocity. 

The second main observational effect of a weak preshock field - an expanded cooling 
zone distance - holds some promise as a diagnostic, and has not been explored completely 
from a theoretical or observational standpoint. In non-magnetic shocks, this distance scales 
approximately as a power law, Vlpi -1 , where Vs is the shock velocity, n 0 is the preshock 
density, and the power x ~ 4.0 — 4.5 (Hartigan et al. 1987). Observationally, one can only 


measure the cooling zone size if there is an emission line that peaks at or near the location 
of the shock, another that peaks within the cooling zone, and this separation is resolvable 
with current instrumentation. Shock velocities in stellar jets are well below the threshold 
where radiation ionizes the preshock gas (~ 100 kms -1 ), so a significant fraction of the H 
that enters the shock is neutral. This neutral H becomes excited by collisions at the shock 
and radiates strongly in Balmer emission lines in a process similar to that which produces 


broad and narrow Balmer filaments in some non-radiative supernova remnants (Raymond 


et al. 2008). Hubble images of stellar jets show that Ha; follows sharp arc-shaped features 


that define the locations of the shocks, while emission lines such as [S II] A6731 follow behind 
these shocks with an observable spatial separation (Heathcote et al. 1996). Recent images 
of HH 2 show enhanced Balmer decrements (Ha/H f3) at the shock fronts, as expected for 


collisional excitation of hydrogen (Raga et al. 2015) 


Cooling zone sizes become more difficult to measure when the shock velocity is high (> 
100 kms -1 ), because the nature of the Balmer emission changes at higher shock velocities. 
Postshock temperatures in high-velocity shocks exceed 10 6 K, and the many different species 
of highly-ionized metals present at these temperatures emit energetic photons capable of 
ionizing H. As a result, most of the Balmer emission no longer occurs from collisional excita¬ 
tions of neutral H immediately behind the front, but instead originates from recombination 
at a secondary peak in the cooling zone where Lyman continuum photons emitted by the 
hot gas near the front are absorbed. In these high-velocity shocks, radiation also affects 
the forbidden line fluxes as electrons ejected by the Lyman continuum photons excite the 
low-energy states of metals by collisions. While full radiative transport is included in the 


models we present in Section 222 radiation plays almost no role in the results of this paper 
because, as noted above, the internal shock velocities within stellar jets derived from both 
emission line ratios and differential proper-motion measurements are typically < 80 kms -1 


Hartigan et al. (1994). 


Our models show that this secondary Ha peak (e.g. Raga & Binette 1991) spreads out 
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over the same temperature range as the forbidden lines of [N If] and [S II], and its integrated 
flux is typically less than that of the integrated flux of the Ha at the front as long as the 
preshock ionization fraction X# < 0.5. In practice, this second peak should not confuse the 
interpretation for most stellar jets because it lies on top of the forbidden line emission and 
does not appear as a distinct sharp feature. Once shock velocities exceed ~ 100 kms -1 , and 
the preshock radiation field becomes intense enough to fully preionize the incoming gas, it 
will become impossible to measure cooling zone distances from Ha. However, in those cases 
one could in principle look for high-ionization lines such as those from N V, O VI, and C IV, 
or even X-ray emission to define the spatial location of the shock front and then measure 
this offset relative to the forbidden-line cooling zone. 

Within the cooling zone there are multiple possibilities for emission-line tracers. For 
example, the lines of [O I] and [S II] both peak around 7000K, while those from [N II] emit 
on average a bit closer to the front, where T ~ 10 4 K. In this paper we define a cooling 
zone distance do as the average distance from the shock of photons from the red [S II] lines 
AA6716+6731. In practice, this is effectively equal to the spatial offset between Ha and 
[S II] because Ha peaks immediately behind the shock. As we show in the next section, this 
measurement of d c, combined with the [S II] A6716/A6731 flux ratio, makes it possible to 
infer Alfvenic Mach numbers in low-velocity radiative shocks. This method is independent 
of abundances and reddening, and uses the brightest emission lines in most jets. The method 
can be refined by including the [N I]/[N II] line ratios integrated over the cooling zone as we 
describe below. 


Shock models in 1-D will not represent HH knots accurately if the cooling distances are 
a substantial fraction of the size of the shock. As we shall see in Sec. [3j cooling distances 
in stellar jets are typically 10 - 100 AU. The width of a typical jet at ~ 0.05 pc from the 
source is > 1000 AU, so most HH shocks are at least an order of magnitude larger than 


the cooling distances. In several cases, spatial offsets between Ha and [S II] exist along the 
entire extent of a well-defined shock front. Such objects are ideal for the analysis described 
here. Shocks that fragment into clumps < 100 AU in size will be problematic to analyze 
with this method. Time-resolved observations of some HH knots have shown significant Ha 
flux variations at the shock in response to variation in the preshock density, but the cooling 
zones are more stable, with flux variations on the order of factors of two over a decade and 
no strong evidence for variations in the line ratios there (Hartigan et al. 2011). Because our 
analysis relies upon integrated fluxes over the cooling zone, with Ha used only to identify 
the location of the shock (Ha intensity is not used), the assumption of steady-state is a good 
one as long as the shocks under consideration have had enough time to develop a full cooling 
zone. This cooling time is typically a decade or two, and only a handful of HH shocks are 
young enough for this to be a concern. 







2.2. Shock Models and Diagnostic Diagrams 


2.2.1. Shock Model Grid 


The shock models are based on the Raymond-Cox 1-D steady-state code (Raymond 


1979). Solar-like logarithmic elemental abundances were used, H: He: C: N: O: Ne: Mg: Si: 
S: Ar: Ca: Fe: Ni = 12.00: 10.93: 8.52: 7.96: 8.82: 7.92: 7.42: 7.52: 7.20: 6.80: 6.30: 7.60: 
6.30. The thermal and density structure of the postshock regions do not depend strongly 
upon the abundances, though ratios of line fluxes from different elements vary in tandem with 
their respective elemental abundances. The code follows the non-equilibrium distribution of 
all ionization states of the 13 elements listed above, and includes all the relevant atomic 
physics such as collisional ionization, photoionization, autoionization, recombination and 
dielectronic recombination, and charge exchange. We have incorporated the latest atomic 
parameters for Fe II given by Bautista et al. (2015) into our calculations. The models 
track the populations of each excited state of all the ions through a matrix of Einstein-A 
values and collision strengths for each transition. For example, populations in 159 levels of 
Fe II are derived in each time step after accounting for over 4000 emission line transitions. 
Radiative losses from both emission-line and continuum processes such as 2-photon radiation 
and bremsstrahlung are included in the energy equation to predict how the fluid variables 
change with time. The code employs several step sizes designed to track any regions of 
abrupt changes in the ionization states of the elements. We chose these as appropriate for 
each of the models we computed. 


The radiative transfer is performed iteratively as follows. At the start of the simulation, 
we use a radiation field from a previous shock model with similar input parameters, and 
track how this spectrum is absorbed downstream in the postshock gas. Meanwhile, all the 
radiative emission lines are stored in an output radiation hie, binned in 1.0 eV intervals 
beginning at 3.6 eV with 150 bins in all. At the end of the simulation the output radiation 
hie is used as the input for the next iteration. Typically this procedure converges within 
two or three iterations to the final model, which predicts huxes for all the emission lines and 
generates temperature, velocity, and density prohles throughout the cooling zone. 


The grid consists of 8470 models. There are 11 choices for the preshock number density 
of nucleons n 0 (cm -3 ), with log(n 0 ) = 2.0 - 4.0 in increments of 0.2. The ten choices for 
preshock ionization fraction X# = n 0 (HI)/n 0 (HI + HII) of H were 0.01, 0.03, 0.05, 0.10, 0.15, 
0.20, 0.30, 0.40, 0.60, and 0.80. Preshock He is neutral in all models. The shock velocity 
Vs had 11 values, starting at 30 kms -1 and ending with 80 kms -1 with increments of 5 
kms -1 . Each set of (Vs, n D , and X#) used seven different preshock magnetic field strengths, 
designed to produce Alfvenic Mach numbers M^ = Vs/Va spaced regularly in log, with 
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log 10 M^ = 1/6, 1/3, 2/3, 1, 4/3, 5/3 and 2 (M A = 1.5, 2.2, 4.6, 10, 22, 46, and 100). The 
Alfven speed is defined by the equation Y\ = B^/(47rn 0 m), where B 0 is the component of 
the preshock magnetic field that lies perpendicular to the shock-normal, n Q is the preshock 
number density of nucleons, and m is the average mass per nucleon. Step sizes are taken 
to be small enough so that the temperature changes by less than 5% between steps when T 
> 10 4 K. Step sizes in the lower-temperature cooling zone (T < 10 4 K) were chosen small 
enough so that the uncertainties in the relevant line ratios introduced by the finite grid size 
were < 1%. Models were terminated when the temperature fell below 2500 K. The preshock 
temperature was 10 4 K in all models. This value is of little consequence in most of the 
models because it represents a small fraction of the kinetic energy of the flow incident upon 
the shock. The sound speed for mostly neutral gas at 10 4 K is only ~ 12 kms -1 . 


2.2.2. Cooling Zone Emission Structure and a New Diagnostic Diagram for Magnetic 

Fields in Shocks 

Fig.0 displays the typical structure present in a shock with a very weak magnetic 
field. The gas in this 50 kms' 1 model enters the shock with a low ionization fraction of 
3%, and reaches a postshock temperature of about 77000 K. Rapid cooling, due in large 
part to collisional excitation of neutral H, quickly drops the temperature to 10 4 K after 
a distance of 0.1 AU, with a more gradual decline that extends over a cooling zone of ~ 
2 AU. The ratio of magnetic to thermal pressure rises as the gas cools, but remains less than 
unity throughout, so the magnetic field has no effect on the gas dynamics in this case. The 
density n rises as the temperature drops to keep approximate pressure balance, producing a 
total compression of about 70 at 2 AU, where the ionization fraction n e /n (inferred from the 
separation of the dashed and dotted curves in the lower panel) is about 10%. The spatial 
distribution of the line emission in the upper panel shows a clear offset between Ha, which 
occurs at the shock, and the forbidden lines. The average distance of the [S 11] emission from 
the front is 0.72 AU. 

A magnetic case is shown in Figs. [2] and [3} In this model, the preshock magnetic 
field is 174 pG, which results in an Alfvenic Mach number of 4.6. The magnetic pressure 
dominates the gas pressure beyond about 0.1 AU, so the density remains nearly constant in 
the cooling zone, with a compression n/n 0 of about 6, and an ionization fraction of ~ 10%. 
This compression is much lower than in the weaker magnetic model in Fig. [l] The field 
reduces the immediate postshock temperature slightly relative to the nonmagnetic case, but 
the main difference between the two cases is that the field extends the cooling zone size by 
a factor of about 20. The average [S If] photon now emits at 19 AU from the shock front. 
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The new models are in agreement with past work (Hartigan et al. 1994) which showed 
that emission line ratios integrated over the entire cooling zone are degenerate with respect to 
field strength in the sense that line ratios from a shock with a small magnetic field resemble 
those from a shock with a stronger field, higher preshock density and higher shock velocity. 
However, when we include the cooling zone distance dc as a constraint, and instead of the 
field strength we use the Alfvenic Mach number as the magnetic observable, we obtain the 
remarkable result shown in Fig. |4j The figure shows that different Alfvenic Mach numbers 
separate in a plot of the [S II] A6716/A6731 line ratio integrated over the cooling distance 
against dc- As both the [S II] ratio and dc are straightforward to measure for any resolved 
cooling zone, it is easy to use this diagram to measure M 4 from observed spectra. 


2.2.3. The Physics Behind the New Diagnostic 

While the results displayed in Fig. [4] can be taken at face value as arising from a re¬ 
liable and well-tested model that includes all the relevant physics of radiative shocks, it is 
a worthwhile exercise to try to understand the main physics that underlies the reason why 
models with different Alfvenic Mach numbers separate in Fig. [4] Let us denote the preshock 
gas with the subscript ‘o’, the gas immediately behind the shock with subscript ‘ 1 ’, and the 
point in the cooling zone where the magnetic and thermal pressures are equal with subscript 
‘2’. The Alfvenic Mach number M 4 is defined in terms of the shock velocity Vs, mass per 
nucleon m, preshock density n D , and preshock magnetic field B 0 as 


M\ = 


4:7T7i 0 mVg 


Bl 


( 1 ) 


Equating the thermal and magnetic pressure at position 2, and using the fact that B ~ n 
throughout the shock we obtain 


Bl 

87 r 


Bln\ 

87ml 


n 2 kT 2 . 


( 2 ) 


If we now make the approximation that the pressure between points 1 and 2 is constant, and 
employ the jump conditions for a strong shock we have 


m =nikTi= ^Yi 

8716 


-n 0 mVl 


(3) 


Combining equations [l] and [3] we get 
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Ma = (4) 

V 3 n 0 

In the magnetic-dominated regime the density remains approximately constant because 
the total pressure ~ B 2 is constant and B ~ n. Hence, equatioii[I]iniphes that the compression 
in the forbidden-line emission zone (relative to the preshock density) is ~ 1.2 M^. For 
example, in Fig. [l] the compression n/n D ~ 6 in the forbidden emission zone, with M 4 ~ 5. 

We now use these relations to help understand the points in Fig. [4} The y-axis of Fig. [4] 
is the flux ratio of [S II] A6T16/A6T31. As is well-known from ISM physics (e.g. 

, there are five low-lying states for any p 3 electronic configu¬ 
ration, and the relative populations of the 2 D 5 / 2)3 / 2 states depend upon the electron density. 
Between the low-density limit (LDL, n e ~ 20 cm -3 ) and the high-density limit (HDL, n e ~ 
2 x 10 4 cm -3 ), Ia 6716 /Ia 673 i scales roughly linearly with log 10 (n e ), with a slope of ~ —0.5. 
Assuming for the moment that cooling zones have the same ionization fractions, the y-axis 
in Fig. [4] becomes 


fe Ferland||2006 Hartigan||2008 ) 


Osterbrock 


y = A6716 ~ constant — 0.5 log 10 n 2 

4A6731 


(5) 


We now consider the cooling distance do Let the volume emission coefficient be e(T) 
(ergcm 3 s -1 ). Balancing the ram energy flux incident upon the shock (ignoring incident 
thermal and magnetic energies, which are small) with the cooling we have 


rd c 


-mn 0 Vg 


e(T(x))n 2 (x)dx ~ en 2 dc 


( 6 ) 


The average emission coefficient e is a complex function of the temperature, and hence 
Vs, because different ionization states of different elements dominate the cooling at distinct 
zones throughout the shock. I 11 the non-magnetic case when shock velocities exceed ~ 
100 kms -1 , numerical models show that d^ ~ n -1 VL where q ~ 4.2 — 4.6 (McKee & 


Hollenbach 1980 Hartigan et al. 1987). In these shocks, the postshock temperature exceeds 


10° K, a regime where equilibrium cooling curves show that e declines with T (Raymond & 


Smith 

1977. 

Schure et al. 

2009) 


material, so a reasonable estimate for the average density in the volume is one that is 
proportional to the postshock density immediately behind shock, which in turn is a factor 
of four higher than the preshock density. Hence, we can estimate n to be proportional to 
For e rjp/2 

r\j V p s , equation [6] implies dc ~ n D 1 r . The numerical simulations 


n, 




















12 


imply an average p ~ —1.5 for T > 10 5 K, consistent with the characteristic slope for these 
temperatures in the equilibrium cooling curves. 

The low-velocity magnetized shocks we are considering in this paper behave differently 
from their high-velocity non-magnetic counterparts for two reasons. First, the shock veloc¬ 
ities are not high enough to fully preionize H. As a result, gas in low-velocity shocks cools 
rapidly through collisional excitation of H immediately behind the shock (Fig. [I]), so that 
most of the gas in the cooling zone has a temperature significantly lower than that of the 
immediate postshock gas. Second, as shown in Figs. [2] and [3j once the magnetic pressure 
exceeds the thermal pressure, the gas maintains nearly a constant density while it cools. 
Hence, high-temperature material near the shock front occupies an even smaller fraction 
of the total cooling volume in the magnetized case (Fig. [2| than it does for non-magnetic 
low-velocity shocks (Fig. [TJ. 

The above discussion implies that for our suite of models, the average density n in the 
cooling zone is approximated well by n 2 , not by a constant times n 0 as it would be in the 
non-magnetic case. Taking n ~ n 2 and solving for dc in equation [6] we get 

dc ~ f(y s )M^W (7) 

where /(Vs) = V|/e is a function of the shock velocity V 5 , and we have used equation [ 4 ] to 
relate n 0 /n 2 to M^. As above, taking e ~ V p s and /(Vs) ~ Vg , where p is an approximate 
power law index gives 


y = constant + 0.5 log 10 dc — 0.5(3 — p) log 10 Vg + 0.5 log 10 A /4 (8) 

Hence, models with fixed Vs and fixed but different n 0 and B 0 follow an approximately 
linear curve in Fig. [fj with the higher-density models having a shorter cooling zones (dc ~ 
f(Vs)M^ 2 n' 1 from equations [ 4 ] and [7]) and correspondingly lower Ia 6716 /Ia 6731 ratios. Alter¬ 
natively, if we fix Vs and dc allow both n 0 and to vary, the positive sign in front of 
the log 10 Ma term in equation [8] implies that models with larger trend upward in Fig. [4j 
This is the physical reason for the color offsets in that figure. In hindsight, that the models 
in Fig. [4] separate when we sort by the dimensionless parameter of Alfvenic Mach number 
rather than by the magnetic field strength B makes sense physically because determines 
how fast the flow is relative to the magnetic signal speed, and is the factor that determines 
the compression behind the shock. 

The points in Fig. [4] also shift as the shock velocity Vg and the preshock ionization 
fraction change, but these effects are significantly smaller than the magnetic effect and 
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effectively introduce some scatter or ‘width’ to each color band. The shock velocity term in 
Eqn. [8] has an opposite sign to that of the Ma term, and so affects the points in Fig. [4] in the 
opposite sense, in that higher values of Vs move points lower in y. For example, the groups 
of 5 - 10 red points that form ‘strings’ that tend down and to the right in the figure are 
models that have the same preshock ionization fraction, preshock density and Alfvenic Mach 
number, but increasing shock velocities and correspondingly higher magnetic held strengths 
as points trend to the lower right. Gaps between these strings of red and black points simply 
result from the finite number of models run. 

The magnetic offsets in Fig [4] dominate over effects caused by different Vs because Ma 
varies by a factor of 67 in the grid, while V changes by less than a factor of 3. The models 
indicate that p ~ 2 over the range of shock velocities in the grid, so the coefficient in front 
of the logio Vs term in equation [8] is essentially the same as the coefficient in front of the 
logio Ma term. Returning to the preshock ionization fraction, we find that it also adds a 
small amount of scatter to Fig. [4j Higher preshock ionization fractions imply higher electron 
densities and more cooling, which leads to points moving down and to the left in Fig. [4j 
The effective scatter introduced is ~ 0.05 in Ia 6716 /Ia 673 i> compared with a scatter of ~ 0.08 
caused by varying Vs- The scatter is somewhat higher at lower M^. 


2.2-4- Improving the Diagnostic Diagram 

Hundreds of possible choices exist for the y-axis in Fig. |4j, using line ratios or some 
combination of line ratios. The [S II] A6716/A6731 ratio has several advantages in that it 
is independent of the reddening and elemental abundances, and is one of the brightest lines 
in an optical spectrum. However, as described above, degeneracies involving Vs and the 
preshock ionization fraction add some scatter to the plot. The [S II] ratio also compresses 
the points near the HDL and LDL, limiting broader applicability. 

We can significantly improve upon uncertainties in estimates of Ma by adding additional 
line ratio measurements. Ideally, each color band in Fig [4] should be made as narrow as 
possible. A new line ratio is useful in reducing the scatter in Ma in Fig. [4] if there is a trend 
in a plot of [S II] A6716/A6731 ratio vs. the new ratio for a fixed range of d<s and a given 
Ma- Mathematically, it is the equivalent of searching for a nonzero partial derivative for a 
new line ratio with respect to [S II] A6716/A6731, keeping dc and Ma fixed. 

Keeping the constraint of bright emission lines that don’t depend upon reddening, and 
recognizing that we need a line ratio sensitive to the ionization, the [N II] /Ha ratio integrated 
over the cooling zone emerges as a natural secondary diagnostic. With [S II] already a 
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diagnostic of density and dc, it is probably not too surprising that an important secondary 
indicator of could come from [N If] A6583. This line traces the ionization fraction of 
N because the N II/N 1 ratio is tied to that of H II/H 1 through a strong charge exchange 
coefficient in the plasma. Another ratio that showed promise was [S II] AA6716+6731/Ha. 


However, the problem with [N II]/Hct, [S II] AA6716+6731/Ha, and similar ratios that 
involve different elements is that they depend linearly upon the elemental abundance ratios. 
These ratios are typically unknown in star-forming regions, and there is even evidence that 
gas-phase abundances of some elements vary along stellar jets as dust grains are destroyed 
(Giannini et al. 2004). Abundance ratios that use different ionization stages of the same 


element avoid this problem, but unfortunately there are no bright [S I] or [S III] lines to 
complement [S II] in most stellar jets, and one cannot easily use [O I]/[0 II] either because 
the brightest [O II] lines (AA3727+3729) are in the blue and subject to large and poorly- 
known reddening corrections, and the red [O II] doublet lines at ~ A7325 are typically faint 
and often blended with a [Ca II] line (e.g. Brugel et al. 1981). One promising line ratio 
remains: [N I] AA5198+5200 / [N II]AA6548+6584, which we henceforth designate simply as 
‘N I/N IT. These [N I] and [N II] lines are typically bright in HH jets, differential reddening 
corrections are small, and abundance uncertainties are not an issue. 


We can visualize the the constraints of additional line ratios geometrically by plotting 
them along a third dimension. In this case we rotate the cube of (log 10 dc, [S II] A6716/A6731, 
and N I/N II) values in such a way as to achieve maximum separation between surfaces of 
constant M^. Fig. [5] shows that models with constant define a series of nested curved 
sheets that become degenerate only in the region where dc < 10 AU, a regime unresolvablc 
spatially with existing instrumentation in any case. There are several orientations that 
separate the different M^, but a particularly convenient one is when the line of sight is in 
the plane defined by (log 10 dc, N I/N II) angled about 18 degrees from the N I/N II axis. 
This projection, which has [S II] A6716/A6731 along the x-axis, and a linear combination of 
log 10 dc and N I/N II along the y-axis, is shown in Fig. [ 6 ] The width of the bands defining 
the red and orange points (M^ = 2.2 and 4.6, respectively) narrows by about a factor of two 
compared with Fig. [4] Most of the remaining scatter in each color arises because surfaces of 
fixed M 4 are not planar sheets, but curve slightly at the ends (Fig. [5]). 

Fig. [ 6 ] provides a quick way to plot a given observation on a graph and determine M^, 
but given a table of models one can also simply collect all the models that are consistent 
with the data. As an example, Fig. [7] analyzes a hypothetical object with three observables, 
dc, [S II] A6716/A6731, and N I/N II. A total of 1131 out of the 8470 models are consistent 
with the constraint [S II] A6716/A6731 = 0.8 ± 0.1. That number drops to 176 models once 
the constraint of [N I]/[N II] = 0.18 ± 0.04 is included, but the consistent models remain 
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spread out over a range of M^. The cooling-zone observation of log 10 d c = 1-1 ± 0.1 on its 
own limits the number of applicable models to 525, but these are again spread out over a 
range of M^. ft is only when the measurement of d c is combined with the [S 11] ratio (163 
models) that a trend in emerges, and the trend becomes a narrow range of values (21 
models) once both [N II]/[N I] and the [S 11] ratio are combined with dc- One drawback 
with such histograms is that the number of models consistent with a given set of constraints 
depends upon the grid used to generate the models. Nevertheless, the peak value of the 
histogram always corresponds to the indicated by the color plot in Fig. [6] 


3. Application to Stellar Jets 


While several examples of resolved cooling-zone distances have been observed with HST, 
images of the ratio of [S 11] A6716/A6731 are rare, and none seem to be available for the 
objects like HH 111 that have relatively simple, well-defined shocks that move in the plane 
of the sky (which minimizes projection effects). Moreover, [N If] images are rare, and [N 1] 
AA5198+5201 images even more rare, despite the fact that both sets of lines are easily 
observed in long-slit spectra. 


For the purposes of this paper, we choose the isolated bow shock ‘K’ in HH 111 for 
further analysis (Reipurth et al. 1997). This object has a resolved cooling zone of log 10 dc = 
1.86 ± 0.14 from archival HST images. Unfortunately, we must attempt to recover line ratios 
from a mixture of ground-based and space-based spectra. Long-slit STIS spectra smoothed 
over 0.4" published by Gonzalez-Gomez & Raga (2003) indicate an electron density of 9000 
cm -3 , but the S/N of the line ratio from these data is very low. From the ground, Morse 


et al. (1993) inferred N e ~ 600 cm” 3 near knot K, which corresponds to [S If] A6716/A6731 
~ 1.0. However, the knot is blended to some degree with knot L, where [S If] A6716/A6731 
~ 1.20. We adopt 1.0 for [S 11] A6716/A6731 in knot K, but the uncertainty is rather large, ~ 
0.20 for this ratio. Taking the Morse et al. (1993) values integrated over knots D through J, 
the observed ratio of [N I] AA5198+5200 / [N II] A6583 is 1.30, and the reddening-corrected 
ratio is 2.49. From the ratio of the A-values, the [N II] A6548 flux should equal 1/3 that 
of [N II] A6583. Hence, log 10 (N I/N II) = —0.01 without reddening corrections, and 0.27 
with reddening included. If we use the observations for knot L, then the reddening-corrected 
value of log 10 (N I/N II) = 0.15. The [N I]/H f3 flux ratio for knot L published by Morse et al. 


(1993) is in agreement with that of Noriega-Crespo et al. (1993), though we cannot use the 
latter reference to estimate [N II] fluxes because those spectra are only at blue wavelengths. 
We adopt log 10 (N I/N II) = 0.27 ± 0.25 for knot K. 

The allowed ranges for the above values define the boxed regions in Figs.[4]an d[6} Despite 
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the large uncertainties (which can be reduced significantly with dedicated observations), it 
is clear from both figures that the Alfvenic Mach number will lie along the lower boundary 
of the models, with ~ 1.5 — 2.2. A total of 58 models yield values for d c, N I/N II 
and [S II] A6716/A6731 consistent with the observational constraints, 25 with = 2.2 
and 33 with = 1.5. Hence, the magnetic field is strong enough to play an important 
role in the dynamics of knot K in this jet. Shock velocities, log(n 0 ), B 0 , and X# are all 
poorly-constrained given the relatively large observational errorbars. 


4. Summary 

The main point of this paper is to show that a key observable - the cooling zone distance 
defined by the separation of Ha from forbidden line emission - makes it possible to infer 
Alfvenic Mach numbers (defined as the shock velocity divided by the preshock Alfven speed) 
in low-velocity shock fronts that cool radiatively. The physics behind why the method works 
is straightforward and the effect is unavoidable: even a weak preshock magnetic field is 
compressed at the shock front and continues to rise behind the front as the gas cools and the 
density increases. Once the magnetic pressure becomes comparable to the thermal pressure, 
the cooling zone lengthens. 

Observational requirements include being able to measure the spatial offset between the 
shock and its cooling zone, for example, from a spatial offset between Ha and [S II], and a 
measure of the red [S II] emission line ratio integrated over the forbidden-line emitting zone. 
If the fluxes of [N I] AA5200 doublet and [N II] A6583 line are available, the models yield 
increasingly more precise results. The ideal dataset has not been acquired for this type of 
analysis, but initial application of the procedure to a bow shock in the HH 111 stellar jet 
indicates a rather low Alfvenic Mach number of ~ 2. This Alfvenic Mach number is relative 
to the shock velocities in the jet. The Alfvenic Mach number relative to the bulk flow speed 
will be about an order of magnitude higher. 


This research is supported by the Department of Energy National Laser Users Facility 
(NLUF) grant de-na0002037, 
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Fig. 1.— Structure of the cooling zone behind a typical low-velocity radiative shock. Top: 
Spatial distribution of radiation from several emission lines, normalized to their peak values 
and plotted as a function of distance from the shock. Bottom: Physical conditions in the 
postshock region. Temperature is read from the scale at left, and the ratios n/n 0 , n e /n 0 , 
and P_e/P ffas refer to the scale at the right. Here, n D is the preshock number density of 
nucleons, n and n e are the number densities of nucleons and electrons, respectively, in the 
postshock region, and P b / P ff as is the ratio of the magnetic pressure to gas pressure. The 
model parameters of Vs, X D , B 0 , n 0 , and M 4 refer to the shock velocity, preshock ionization 
fraction, preshock magnetic held strength, preshock density, and incident Alfvcnic Mach 
number, respectively. The variable dc is the [S II]-weighted cooling-zone distance, defined 
in the text. This model is nonmagnetic in the sense that magnetic pressure is less than the 
gas pressure everywhere. 
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Fig. 2.— Same as Fig. [I] but for a stronger preshock magnetic field. The magnetic held 
pressure expands the cooling zone distance relative to the nonmagnetic case. 
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Fig. 3.— Same as Fig. [2] with an expanded scale on the x-axis. 
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Fig. 4.— Different Alfvenic Mach numbers plotted as a function of two observables, the 
[S II] A6T16/A6T31 emission line ratio and the cooling distance, easily measured from the 
spatial offset between Her and [S II] emission. The separation of Alfvenic Mach numbers 
in this diagram provides a means to estimate the importance of magnetic fields in the flow 
dynamics. The boxed area represents current observational uncertainties for knot K in 
HH 111. 
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Fig. 5.— Same as Fig. [4] but adding in a third dimension N I/N II, defined as [N I] 
AA5199+5201/[N II] AA6548+6583. Models with fixed trace nearly parallel sheets. 
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Fig. 6.— Same as Fig. [|but including a correction factor to dc that is based on the observed 
N I/N II ratio. This graph is the projection of Fig. [5] along a line of sight oriented 18° from 
the N I/N II axis. The boxed region denotes the possible range of values related to knot K 
in HH 111. 
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Fig. 7.— Histograms of the number of models satisfying the listed constraint for a hypothet¬ 
ical object. Alfvenic Mach numbers require both line-ratio measurements and an estimate 
of the cooling distance dcoon- 
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